{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imaging with a vector-Apodizing Phase Plate coronagraph\n",
    "\n",
    "We will introduce the propagation of polarized light through a vector apodizing phase plate coronagraph (vAPP).\n",
    "\n",
    "We'll start by importing all relevant libraries and setting up our pupil and focal grids. We also import precomputed phase file for the pupil and the vAPP and the pupil the vAPP pattern was calculated for. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hcipy import *\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "# For notebook animations\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pupil_grid = make_pupil_grid(512)\n",
    "focal_grid = make_focal_grid(4, 40)\n",
    "\n",
    "prop = FraunhoferPropagator(pupil_grid, focal_grid)\n",
    "\n",
    "vAPP_phase = Field(read_fits('vAPP_phase.fits.gz').ravel(), pupil_grid)\n",
    "telescope_pupil = Field(read_fits('vAPP_amplitude.fits.gz').ravel(), pupil_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A vector apodizing phase plate can be modelled as a half-wave retarder with a varying fast-axis orientation [1] [2]. A vAPP applies geometric phase to the incoming wavefront. The applied geometric phase, $\\delta (x, y)$,  only depends on the fast-axis orientation, $\\Phi (x, y)$, through:\n",
    "\n",
    "\n",
    "$$\\delta (x, y)= \\pm 2 \\Phi (x, y).$$\n",
    "\n",
    "The sign of the acquired phase depends on the input circular polarization state. To simulate a vAPP, we define it as a phase retarder with a fast-axis orientation that is the phase / 2 and a retardation of $\\pi$.\n",
    "\n",
    "We plot the phase pattern and the fast-axis orientation pattern. \n",
    "\n",
    "[1] Snik, Frans, et al. \"The vector-APP: a broadband apodizing phase plate that yields complementary PSFs.\" Proc. SPIE, Vol. 8450 (2012)\n",
    "\n",
    "[2] Otten, Gilles PPL, et al. \"Performance characterization of a broadband vector Apodizing Phase Plate coronagraph.\" Optics Express 22.24 (2014)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quiver_plot_retarder(fast_axis_orientation, N):  \n",
    "    '''Plot the fast-axis orientation of a retarder\n",
    "\n",
    "    Parameters\n",
    "    ---------\n",
    "    fast_axis_orientation : Field\n",
    "        The input fast_axis_orientation that is plotted\n",
    "    N : Scalar\n",
    "        A scalar that determines the sampling of fast_axis_orientation that is plotted\n",
    "    '''\n",
    "    pupil_grid = fast_axis_orientation.grid\n",
    "    X = pupil_grid.x.reshape(pupil_grid.shape)[::N,::N]\n",
    "    Y = pupil_grid.y.reshape(pupil_grid.shape)[::N,::N]\n",
    "    U = np.cos(fast_axis_orientation.shaped[::N,::N])\n",
    "    V = np.sin(fast_axis_orientation.shaped[::N,::N])\n",
    "    quiveropts = dict(color='black', headlength=0, pivot='middle', scale=30, units='xy', width=.003, headwidth=0)\n",
    "    plt.quiver(X,Y,U,V,**quiveropts )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the properties of the vAPP\n",
    "fast_axis_orientation = vAPP_phase / 2\n",
    "phase_retardation = np.pi\n",
    "circularity = 0\n",
    "\n",
    "# Making the vAPP \n",
    "vAPP_pol = PhaseRetarder(phase_retardation, fast_axis_orientation, circularity)\n",
    "\n",
    "# Setting up the sampling of the quiver plot\n",
    "N = 16\n",
    "\n",
    "# Plotting the phase pattern and fast-axis orientation pattern\n",
    "plt.figure()\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "imshow_field(fast_axis_orientation, mask=telescope_pupil, cmap='RdBu')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "imshow_field(telescope_pupil, cmap='gray')\n",
    "quiver_plot_retarder(fast_axis_orientation, N)\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vAPP phase pattern has a phase ramp ( = polarization grating) to separate the two circular polarization states when imaged on a detector. The sampling of the fast-axis orientation is not good enough to fully capture the phase pattern. Therefore, we subtract the grating pattern from the phase pattern. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Removing the grating pattern\n",
    "fast_axis_orientation_new = (vAPP_phase - pupil_grid.y * 20 * 2 * np.pi) / 2\n",
    "\n",
    "# Changing the values to make the pattern aesthetically pleasing\n",
    "fast_axis_orientation_new = (fast_axis_orientation_new - np.pi / 2) % np.pi\n",
    "\n",
    "# Setting up the sampling of the quiver plot\n",
    "N = 16\n",
    "\n",
    "# Plotting the phase pattern and fast-axis orientation pattern\n",
    "plt.figure()\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "imshow_field(fast_axis_orientation_new, mask=telescope_pupil, cmap='RdBu')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "imshow_field(telescope_pupil, cmap='gray')\n",
    "quiver_plot_retarder(fast_axis_orientation_new, N)\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate the vAPP by propagating a left-circularly polarized wavefront through the vAPP. The intensity image is simulated by propagating to the focal plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavefront = Wavefront(telescope_pupil, 1, input_stokes_vector=(1, 0, 0, -1))\n",
    "wavefront_out = vAPP_pol.forward(wavefront)\n",
    "\n",
    "vAPP_PSF = prop(wavefront_out).intensity\n",
    "\n",
    "imshow_field(np.log10(vAPP_PSF / vAPP_PSF.max()), vmin=-5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vAPP PSF has a dark zone from 2-15 $\\lambda/D$ and the coronagraphic PSF is offset from the center because of the grating that is added to the phase. For the right-circular polarization state, the applied phase is opposite, and the PSF is different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavefront = Wavefront(telescope_pupil, 1, input_stokes_vector=(1, 0, 0, 1))\n",
    "wavefront_out = vAPP_pol.forward(wavefront)\n",
    "\n",
    "vAPP_PSF = prop(wavefront_out).I\n",
    "\n",
    "imshow_field(np.log10(vAPP_PSF / vAPP_PSF.max()), vmin = -5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate the observation of a star, we use unpolarized light. Unpolarized light contains equal amount of left-circular polarization and right-circular polarization. Therefore, the stellar PSF is a combination of the two previously simulated PSFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavefront = Wavefront(telescope_pupil, 1, input_stokes_vector=(1, 0, 0, 0))\n",
    "wavefront_out = vAPP_pol.forward(wavefront)\n",
    "\n",
    "vAPP_PSF = (prop(wavefront_out).I)\n",
    "\n",
    "imshow_field(np.log10(vAPP_PSF / vAPP_PSF.max()), vmin = -5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In an ideal case is the vAPP described by a patterned half-wave retarder. However, the retardation is not always half wave. This does not change the coronagraphic PSF, but rather changes the amount of light that is diffracted as the coronagraphic PSFs. The fraction that is not diffracted is called polarization leakage and is imaged in the center of the detector. Here we simulate the image plane of a vAPP for different retardations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "retardation = np.linspace(0, 2 * np.pi, 37, endpoint=True)\n",
    "wavefront = Wavefront(telescope_pupil, 1, input_stokes_vector=(1, 0, 0, 0))\n",
    "vAPP_pol = PhaseRetarder(retardation[0], fast_axis_orientation, circularity)\n",
    "\n",
    "fig = plt.figure()\n",
    "im = imshow_field(focal_grid.x, focal_grid, vmin=-5, vmax=0)\n",
    "title = plt.title('', fontsize='xx-large')\n",
    "plt.close()\n",
    "\n",
    "\n",
    "def animate(i):\n",
    "    vAPP_pol.phase_retardation = retardation[i]\n",
    "    wavefront_out = vAPP_pol.forward(wavefront)\n",
    "\n",
    "    vAPP_PSF = Field((prop(wavefront_out).I),focal_grid)\n",
    "    vAPP_PSF /= vAPP_PSF.max()\n",
    "\n",
    "    x, y = focal_grid.coords.separated_coords\n",
    "    im.set_data(x, y, np.log10(vAPP_PSF.shaped))\n",
    "    title.set_text('Retardance = {0} degrees'.format(round(np.degrees(retardation[i]),0)))\n",
    "    return im, title\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, range(len(retardation)), blit=True)\n",
    "HTML(ani.to_html5_video())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vAPP applies gemeometric phase that is inherently independent of wavelength. The shape of the coronagraphic PSFs does not change as function of wavelength, except for wavelength scaling. In addition, the retardance changes as wavelength, causing the intensity of the polarization leakage to change accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelengths = np.linspace(0.9, 1.1, 21)\n",
    "\n",
    "retardation = lambda wavelength: np.pi / wavelength\n",
    "\n",
    "fig = plt.figure()\n",
    "im = imshow_field(focal_grid.x, focal_grid, vmin=-5, vmax=0)\n",
    "title = plt.title('',fontsize = 'xx-large')\n",
    "plt.close()\n",
    "\n",
    "vAPP_pol = PhaseRetarder(retardation, fast_axis_orientation, circularity)\n",
    "\n",
    "def animate(i):\n",
    "    wavefront = Wavefront(telescope_pupil, wavelength=wavelengths[i], input_stokes_vector=(1, 0, 0, 0))   \n",
    "    wavefront_out = vAPP_pol.forward(wavefront)\n",
    "\n",
    "    vAPP_PSF = Field((prop(wavefront_out).I),focal_grid)\n",
    "    vAPP_PSF /= vAPP_PSF.max()\n",
    "\n",
    "    x, y = focal_grid.coords.separated_coords\n",
    "    im.set_data(x, y, np.log10(vAPP_PSF.shaped))\n",
    "    title.set_text('Normalized wavelength = %.2f' % wavelengths[i])\n",
    "    return im, title\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, range(len(wavelengths)), blit=True)\n",
    "HTML(ani.to_html5_video())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "level": "beginner",
  "thumbnail_figure_index": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
